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SUMMARY 

We introduce a novel nonlinear seismic imaging method based 
on model order reduction. The reduced order model (ROM) is 
an orthogonal projection of the wave equation propagator op¬ 
erator on the subspace of the snapshots of the solutions of the 
wave equation. It can be computed entirely from the knowl¬ 
edge of the measured time domain seismic data. The image is 
a backprojection of the ROM using the subspace basis for the 
known smooth kinematic velocity model. The implicit orthog- 
onalization of solution snapshots is a nonlinear procedure that 
differentiates our approach from the conventional linear meth¬ 
ods (Kirchhoff, RTM). It allows for the removal of multiple 
reflection artifacts. It also enables us to estimate the magni¬ 
tude of the reflectors similarly to the true amplitude migration 
algorithms. 

INTRODUCTION 

To simplify the exposition we make a number of assumptions 
on our model. 

First, we consider an acoustic wave equation 

u tt = Au , fe[o,r], (l) 

for the pressure u. We treat the spatial operator A = C 2 A as a 
matrix A G R NxN , a discretization on some fine grid. The diag¬ 
onal matrix C = diag(c) contains the acoustic velocity values 
c G R n at the N nodes of the fine grid, while A discretizes the 
Laplacian on that grid. 

Second, we assume that the sources can be modeled as an ini¬ 
tial condition 

u(0) = S, 11,(0) = 0, (2) 

which can be easily achieved by considering the even part of 
the time domain solution of equation [I] and thus the even part 
of the data. The matrix S G R Nxp contains all p sources and is 
localized near the surface. 

Under the initial conditions of equation[2]the solution to equa¬ 
tion [l]takes the form 

u (t) = cos(rV—A)S, (3) 

where the cosine and the square root are understood as matrix 
functions. Note that the matrix function of time u (t) G R Nxp 
consists of p columns that contain the solutions corresponding 
to each source in S. 

Third, we suppose that the source matrix S admits a represen¬ 
tation 

S = tf 2 (A)CE, (4) 

where E G R Nxp are p point sources supported on the surface 
and q 2 ((o) is the Fourier transform of the source wavelet. Here 
we consider a Gaussian source wavelet 

q 2 (A) = e aA . 


For small <7 all quantities in equation]?] are localized near the 
surface. We assume that the velocity near the surface is known, 
thus we know S. 

Fourth, we assume that the sources and receivers are collo¬ 
cated. This assumption makes the construction of the reduced 
order model in the next section a lot more straightforward. 
However, our approach can be generalized to the much more 
realistic case of noncollocated sources and receivers. The par¬ 
ticular form of the receiver matrix R G R Nxp that we use is 

R = C -1 E. (6) 

Finally, under all the assumptions above we can write a data 
model 

F(f) =R r cos(;\/ = A)S, (7) 

which is a px p matrix function of time. 

The problem of seismic imaging that we solve here is to find 
an estimate of the unknown velocity c from the knowledge of 
the time domain data F(f) for t G [0,7] and a smooth kinematic 
velocity model Cq. 

TIME-DOMAIN INTERPOLATORY REDUCED ORDER 
MODEL 

At the core of our approach is the construction of the ROM that 
interpolates the measured seismic data. The use of model or¬ 
der reduction techniques in inversion was proposed in |Druskin| 
|et al.| ( |2013| >; |Borcea et al.| ( |2014| ) for parabolic (controlled 
source electromagnetic method, CSEM) inverse problems. Un¬ 
like the diffusive parabolic case, where the authors employed 
frequency (Laplace) domain interpolation, the appropriate set¬ 
ting for the wave equation inversion is the time domain. 

The particular form of the source and receiver matrices in equa¬ 
tions |?J|6] allows us to rewrite the data model from equation [7] 
in the completely symmetric form 

F(f)=B r cos()\/-T)it (8) 

where the symmetrized spatial operator is A = CAC and the 
source/receiver matrix is given by B = g(A)E. Here we used 
the fact that analytic matrix functions commute with similarity 
transforms and also that the symmetric operator A is similar to 
the original A with a similarity transform C. 

In practice, the time domain data is measured at discrete time 
instants that we denote by = kx with k = 0,1,..., 2n — 1, 
where X is the sampling interval and U«-l =7 is the terminal 
time. 

The discrete data samples = F(^) admit a representation 

= F (kx) = B t cos ^k arccos cos x— A^ B = B r 7^(P)B, 

(9) 


(5) 
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where T\ are Chebyshev polynomials of the first kind and P is 
the propagator given by 


P = cos 



( 10 ) 


Applying equation [T7] twice we can also obtain 
(U r PU), 7 = ufPu; = 

= 4 ( F <+J+1 + F |i-;+l| + F |i+j- 1| + F |i-;-l|)- 


( 19 ) 


We wish to construct a reduced order model of size np CiV 
that matches the 2 n measured data samples exactly 

Fj; = B r 7^(P)B = B r 7^(P)B, k = 0,1,.. ,,2n- 1, (11) 

where P G R n P xn P and B G R ni?Xi? . Since we are^solving the 
inverse problem of seismic imaging the ROM (P, B) should be 
computable from the knowledge of the sampled data only. 

It appears that the solution to the data interpolation problem of 
equation [IT] can be found in the projection form 

P = V r PV, B = V r B, (12) 

where the columns of V € M Nxnp constitute an orthonormal 
basis for the subspace spanned by the discrete time snapshots 
of solutions 


u k = u(t/t) = cos (toy --Aj B = 7j.(P)B € R Nx P. (13) 

If we introduce the matrix of solution snapshots 

U = [u 0 ,ui,• • ■ ,u„_i] G R Nxnp , (14) 

then V is defined simply by 

colspan V = colspan U, V r V = I. (15) 

Note that the above definition is not unique as V is defined up 
to an orthonormal change of variables in the projection sub¬ 
space. For the purpose of imaging some choices of V are better 
than others. 

The orthogonalization of snapshots U must respect the causal¬ 
ity and the propagating nature of the time domain solutions of 
the wave equation. Thus, each snapshot should be orthogonal- 
ized only against the previous ones. In linear algebra this is 
known as Gram-Schmidt procedure or the QR decomposition. 
Since our snapshots are matrices with p columns correspond¬ 
ing to all sources/receivers, we need a block version of QR 
decomposition 

U = VL r , (16) 

where \J G M npxnp is block upper triangular with blocks of 
size p. 

Obviously, we cannot simply use equation |T6] since the snap¬ 
shots U are unknown to us. However, from the data we can 
obtain the inner products between the snapshots. A basic mul¬ 
tiplication property of Chebyshev polynomials 

Ti(x)Tj(x) = l ( T i+j (x ) + 2| M (*)) (17) 

combined with equation [13] immediately implies that 

(U r U)jy = uf Uj = I(F i47 - + F| M ). ( 18 ) 


The knowledge of Gram matrix U r U from equation 


18 


allows 


us to compute the block lower triangular factor L in equation 
p~6]via a block Cholesky decomposition 


U r U = LL r 


( 20 ) 


Once the Cholesky factor is known we use equation [T9] to ob¬ 
tain the final expression for the ROM 

P = L _1 (U 7 PU)L _r , (21) 

entirely from the sampled data F^. 


BACKPROJECTION IMAGING 

After the reduced order model of equation [2l]is obtained from 
the measured data we need to extract from it the information 
about the velocity c. The first step is to go from the ROM 
for the propagator P to the reduced model for the symmetrized 
spatial operator A = CAC by approximately inverting equation 
[Housing the first two terms in the Taylor’s expansion 

A= 4(P-I) «V r AV. (22) 

T 2 

There are multiple ways Jo obtain the estimate of the velocity 
from the knowledge of A. One may employ optimization to 
solve for c by minimizing the ROM misfit. Such procedure 
is superior to the conventional full waveform inversion (FWI) 
which minimizes the data misfit. However, in this work we are 
interested in a non-iterative imaging algorithm that assumes 
the knowledge of a smooth kinematic background model de¬ 
noted by Co- 

If the projection subspace colspan V is sufficiently rich, then 
the backprojection must be a good approximation of the spatial 
operator 

A^VAV 7- . (23) 

However, we have no direct access to the orthonormal basis 
V. We approximate it with a known basis Vq for the smooth 
kinematic velocity model Co- 

In order to get an imaging formula we also notice that the di¬ 
agonal of A is proportional to the square of the velocity 

c 2 oc diag(A) = diag(CAC), (24) 

where the square c 2 is understood componentwise. Similarly, 
for the difference between the unknown velocity and the kine¬ 
matic model we can write 

<5c 2 = c 2 — Cq diag(A — Aq). (25) 

Replacing the symmetrized operators A and Aq in equation 
[25] with their backprojection approximations from equation [23] 
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and also using the approximation V « V 0 we arrive at the for¬ 
mula 

Sc 2 ocdiag(Vo(A-Ao)v£), (26) 

where A is computed from the measured data and Aq,\q are 
easily found since Cq is known. 

Observe that unlike the conventional imaging approaches (re¬ 
verse time migration, Kirchhoff migration) the formula in equa¬ 
tion [26] is nonlinear in the measured data. This is due to the 
nonlinearity of the block Cholesky decomposition in equation 
[20] and inversion of the Cholesky factor L in equation |2l] The 
nonlinearity that amounts to the implicit orthogonalization of 
solution snapshots U allows for a better quality image. In par¬ 
ticular, the orthogonalization removes the multiple reflection 
artifacts which are otherwise very difficult to deal with using 
conventional linear migration algorithms. This is illustrated in 
Figure [I] where we show a simple synthetic model with two 
layers. The backprojection images the layers correctly while 
suppressing the multiple reflection artifacts that are present in 
the RTM image as ghost layers below the actual ones. 
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Figure 1: Removal of multiples: (a) true velocity c; (b) back- 
projection image c*; (c) RTM image computed as a gradient 
of conventional FWI. Distances are in km , velocities in km/s, 
p =12 sources/receivers are black x. 


Equation[26]can be used in a number of ways to obtain the im¬ 
age. The ambiguity comes from the choice of the proportion¬ 
ality factor. Here we choose it to be the background velocity, 
which leads us to a multiplicative imaging formula 


(27) 


where c* « c is the image, a is a scalar step length and all al¬ 
gebraic operations on the right hand side are performed com¬ 
ponentwise. With the imaging formula from equation [27] at 
hand we can summarize our seismic imaging method in the 
following algorithm. 


Algorithm 1 (Nonlinear ROM backprojection imaging) 


1 . Choose the sampling time interval z and measure the 
discrete time samples of the seismic data with equation 
[7f F k = F(zk)fork = 1,2,...,2/i-L 

2. Using equation\T8\compute the snapshot Gram matrix 
U T U from the data and perform its block Cholesky de¬ 
composition as in equation\20\ 

3. Using equation 19 compute the matrix U r PU from the 


data and use Cholesky factor from step 2 to form the 
ROM P with equation 


21 


4. 


Choose a smooth kinematic velocity model Cq and use 
the block QR decomposition of equation [7b] to com¬ 
pute the orthonormal basis Yq; project Pq on Vq using 


equation 


12 


to obtain the kinematic model ROM Pq. 


5. From the propagator ROMs P and Pq obtain the spatial 


operator ROMs Aq and Aq using equation 


22 


6. With the operator ROMs from step 5 and the orthonor- 


26 


mal basis from step 4 compute 8c 2 from equation 
and use the imaging formula of equation |27| to form 
the final image c*. 


NUMERICAL EXAMPLE: MARMOUSI MODEL 


We evaluate the performance of our method on the synthetic 
data computed for the Marmousi model by Bourgeois et al. 
fl99l) . 


x=6.90 km 



Figure 2: Vertical well log for the Marmousi model at offset 
x = 6.9 km. Depth (horizontal axis) is in km, velocity (vertical 
axis) is in km/s. True c is black o, smooth kinematic velocity 
model Co is red x and the image c* is green □. 

The model is on a 15 m grid with N = 900 x 180 = 162,000 
nodes. The choice of the time interval z is very important for 
our method. To make the orthogonalization procedure well 
conditioned, it should be chosen at a Nyquist rate for the given 
source wavelet. Here we use z = 33.5 ms {n = 35) which cor¬ 
responds to the frequency of about ft) = 15 Hz. The data is 
measured for p = 90 sources/receivers spaced uniformly ev¬ 
ery 150 m. The kinematic model Co is obtained by convolving 
the true velocity c with a Gaussian kernel of width 465 m and 
height of 315 m. 


c* = cq \J 1 + a<5c 2 , 
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Figure 3: Seismic image for a Marmousi model: (a) True velocity c; (b) Smooth kinematic velocity model Co; (c) Nonlinear ROM 
backprojection image c*; (d) Difference c* — cq. All distances are in km , velocities in km/s. The sources/receivers are black x. 


In Figure [5] we show the image c* and the difference c* — Co 
along with the smooth kinematic and true Marmousi models. 
We observe very good recovery of all the model’s features 
down to 2.4 km. The very bottom is not imaged because we 
had to truncate the data sampling at 2n = 70 to avoid the reflec¬ 
tions from the bottom, which at the moment employs reflective 
boundary conditions instead of the PML. 

We also show in Figure [2] a vertical well log. It demonstrates 
that our method performs well not only recovering the loca¬ 
tions of the reflectors but also their strengths. We observe that 
the magnitude of the imaged velocity c* is in good agreement 
with the true model c. In this particular aspect our algorithm 
performs similarly to the true amplitude migration methods. 


CONCLUSIONS AND DISCUSSION 

We introduced a novel nonlinear seismic imaging method based 
on the backprojection of reduced order models computed di¬ 
rectly from the measured time domain data. The results of the 
early numerical experiments with synthetic data for Marmousi 
model show great promise. 

The main issue to be solved to make the method viable for 
the real field data is to remove the assumption that the sources 
and receivers are collocated. This is certainly possible if one 
uses different left (source) and right (receiver) subspaces in the 
projection equation [12] Other possible improvements include 
the implementation of absorbing boundary conditions (PML) 
and more accurate source models. 



































Nonlinear imaging via ROM backprojection 


REFERENCES 

Borcea, L., V. Druskin, A. V. Mamonov, and M. Zaslavsky, 
2014, A model reduction approach to numerical inversion 
for a parabolic partial differential equation: Inverse Prob¬ 
lems, 30 , 125011. 

Bourgeois, A., M. Bourget, P. Lailly, M. Poulet, P. Ricarte, and 
R. Versteeg, 1991, Marmousi, model and data: The Mar- 
mousi Experience, 5-16. 

Druskin, V., V. Simoncini, and M. Zaslavsky, 2013, Solution 
of the time-domain inverse resistivity problem in the model 
reduction framework Part I. One-dimensional problem with 
SISO data: SIAM Journal on Scientific Computing, 35 , 
A1621-A1640. 



